function dq_Fdd = fn_dq(delta_grid,R1,D0,lambda,prob,dist_s,s_min,s_max,dist_d,d_min,d_max,phi,delta_grid_stp)
% Could be made more precise, needs delta large otherwise

% delta_grid = (0:delta_grid_stp:R1*d_max)';
s_hat      = fn_s_hat(R1,lambda,s_min,s_max,phi);

q_F = zeros(length(delta_grid),1);
for j = 1:length(delta_grid)
    delta_i = delta_grid(j);
    
    s_star  = fn_s_star(delta_i,R1,D0,lambda,s_min,s_max,dist_d,d_min,d_max,phi);
    q_F(j)   = dist_s.cdf(s_hat*ones(size(s_star))) + prob*(dist_s.cdf(s_star)-dist_s.cdf(s_hat));
    
end

dq_Fdd  = gradient(q_F,delta_grid_stp); % better than diff()

% dq_grid = delta_grid;

% figure; plot(delta_grid,dq,'Linewidth',2);

% f_dq = @(x) interp1(dq_grid,dq,x,'pchip'); % Alternative
% dq = f_dq(delta_grid)

end
